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ABSTRACT 


i  V* 

-■sqv-jriv-S? 

The  application  of  Richardson  iteration  to  a  symmetric,  but  indefinite 
linear  system  requires  certain  parameters  which  can  be  determined  from  the 
zeros  in  the  error  of  a  certain  best  polynomial  approximant  on  some  set  S 
known  to  contain  the  spectrum  of  the  coefficient  matrix.  It  is  pointed  out 
that  this  error  can  also  be  obtained  as  a  multiple  of  the  extremal  polynomial 
for  the  linear  functional  p  M-  p(0),  and  this  leads  to  an  efficient  Remes 

/  /in 

type  algorithm  for  its  determination.  A  program  incorporating  this  algorithm 
for  the  case  that  S  consists  of  two  intervals  bracketing  zero  is  also  given. 


AMS(MOS)  Subject  Classification  -  65F10,  65D15,  41A10 
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SIGNFICANCE  AND  EXPLANATION 


Sparse  symmetric  (positive)  definite  linear  systems  occur  in  various 
problems,  chiefly  in  the  construction  of  least-squares  approximations  and  in 
the  solution  of  elliptic  partial  differential  equations.  They  are  most 
advantageously  solved  by  the  conjugate  gradient  method.  But  this  method  may 
break  down  if  the  system  fails  to  be  definite,  as  it  might  be  if  constraints 
are  added  to  the  least-squares  problem  or  standard  conditions  on  the  lower 
order  terms  in  the  elliptic  PDE  are  relaxed.  In  such  a  case,  Richardson 
iteration  offers  a  possibly  attractive  alternative.  This  iteration  method 
requires  knowledge  of  some  set  S  certain  to  contain  the  spectrum  of  the 
coefficient  matrix  of  the  linear  system  and  this  can  often  be  supplied.  In 
addition,  it  requires  knowledge  of  the  polynomial  of  smallest  size  on  S 
which  has  a  given  degree  and  takes  the  value  1  at  0  .  If  S  lies  to  one 
side  of  0  ,  then  this  polynomial  is  just  the  Chebyshev  polynomial  for  S 
(normalized  to  have  value  1  at  0)  .  But,  in  the  indefinite  case,  the 
Chebyshev  polynomial  will  not  do. 


EXTREMAL  POLYNOMIALS  WITH  APPLICATION  TO  RICHARDSON  ITERATION 
FOR  INDEFINITE  LINEAR  SYSTEMS 

Carl  de  Boor  and  John  R.  Rice 


1.  The  iteration  problem  Consider  the  linear  system  of  equations  Ax  =  b  and  the 


iteration 

n+1  n  . .  n 

x  -  x  -  a  ( Ax  -b) 
n 

With  «n  :=  x"  -  x  the  error  in  the  n-th  iterate,  we  have 

en  =  (1-a  A)en  1  =  ...  =  H  (1-a.  A)  e°  =  Q  (A)e° 

n-1  .  ,  1-1  n 

1=1 

where  Qn  is  the  polynomial  of  degree  n  which  vanishes  at  l/ctj,  1/a  1  and  is 

1  at  0  .  This  is  Richardson's  (first  order)  iteration,  with  iteration  parameters  .  If 

the  spectrum  of  A  is  known  to  lie  in  some  compact  set  S  ,  then  a  standard  analysis 

suggests  that  one  should  choose  the  parameters  or  so  as  to  minimize 

Q_  -  :=  max  10  (s>| 

n  o  _  n 

ses 

The  resulting  polynomial  Pn  is  then  the  error  in  the  best  Chebyshev  approximation  on  S 

to  1  from  {E.^B.t-*}  .  If  S  is  an  interval  not  containing  the  origin  (hence  A  is 

—  -  -,=1  3 

known  to  be  definite),  then  it  is  well  known  that  a  renormalization  of  Pn  to  make  the 
coefficient  of  tn  equal  to  1  gives  Tn  ,  the  Chebyshev  polynomial  for  the  interval  S  . 
For  this  case,  the  three-term  recurrence  relation  for  the  Chebyshev  polynomials  may  he 
employed  to  build  up  x11  without  the  use  of  the  zeros  of  Pn  .  This  has  the  advantage 
that  the  iterates  x*  so  generated  along  the  way  are  themselves  using  P^  .  This  method  is 
known  as  the  Chebyshev  semi-iterative  method  .  This  variation  requires  some  more  memory  (3 
vectors  rather  than  2  are  used)  and  more  computation  per  step  (since  more  vectors  are 


combined  per  step).  The  conjugate  gradient  method  is  a  further  variation  which,  with  some 
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more  work  per  iteration#  removes  the  dependence  on  the  interval  S  ;  the  mere  knowledge 
that  such  an  interval  exists  suffices  to  show  that  the  error  produced  at  the  n-th  step  is 
of  the  form  Pne°  with  Pn  the  error  in  a  best  approximation  to  1  on  the  spectrum  of 
A  itself. 

The  conjugate  gradient  method  may  run  into  difficulties  when  A  ,  though  symmetric  and 
invertible,  is  not  definite.  See  Paige  and  Saunders  [1975]  for  a  detailed  discussion  and 
some  remedies.  For  this  reason#  Richardson  iteration  becomes  an  attractive  alternative  in 
this  case.  We  now  have  the  spectrum  of  A  contained  in  two  intervals,  with  the  origin 
between  them.  Akhiezer  [1928]  determined  the  Chebyshev  polynomials  for  two  such  intervals 
of  equal  length  and  Lebedev  [1969]  extended  this  technique  to  a  set  S  consisting  of  an 
arbitrary  number  of  intervals  of  equal  length  and  applied  his  result  to  iteration*  See 
Anderssen  and  Golub  [1972J  for  a  translation  of  Lebedev's  paper  and  further  discussions, 
particularly  on  the  important  subject  of  the  order  in  which  best  to  use  the  a^'s  • 

Specifically,  let  S  *  [a,b]  [c,d]  •  For  certain  values  of  a,  b,  c  and  d,  Atlestam 

[1977]  has  obtained  a  representation  of  the  Chebyshev  polynomials  for  S  ,  of  the 
following  form:  Let 

Q(  t )  cos(m[7T+I(  t) )  ) 


with 


and 


l(t) 


t  c  c 

J  (u-r)p(u)du  ,  r  z~  J  up(u)du  /  /  p(u)du 

a  b  b 


p(u)  :=  ( (u-a)(u-b) (u-c) ( d-u) j  ^ 

If  there  arc  integers  m  and  k  so  that  Kb)  =  kir/m  ,  then  Q  is  a  polynomial  of 
degree  m  proportional  to  the  Chebyshev  polynomial  for  S  .  Atlestam  further  shows  that, 
for  any  interval  pair  S,  the  Chebyshev  polynomial  is  of  this  form  but  for  a  slightly 
different  pair  of  intervals,  and  this  difference  goes  to  zero  as  the  degree  goes  to 
infinity.  Her  arguments  can  be  used  to  show  that  in  the  same  way,  for  any  interval  pair 
Z  bracketing  t-he  origin,  the  best  polynomial  Pn  is  of  the  above  form,  but  for  a  slightly 
different  interval  pair.  These  results  can  be  used  to  obtain  sharp  asymptotic  results  on 
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the  degree  of  convergence  of  the  iteration  method,  but  it  is  not  clear  how  useful  the 
representation  is  for  obtaining  the  necessary  iteration  parameters. 

In  the  present  paper,  we  give  what  we  feel  is  a  more  useful  formulation  of  the 
mathematical  problem  underlying  the  determination  of  the  parameters;  well  known  results 
then  establish  existence  and  uniqueness  of  the  solution  of  this  problem  and  characterize 
it.  In  particular,  we  are  led  to  a  Remes  type  algorithm  for  the  determination  of  Pn  , 
whose  zeros  can  then  be  determined  efficiently  by  the  Modified  Regula  Falsi.  For  the 
particular  case  that  S  is  an  interval  pair,  we  present  some  numerical  results  to 
illustrate  the  nature  of  the  parameters  and  the  convergence  rate  of  the  corresponding 
iteration.  In  an  appendix,  we  list  a  Fortran  program  (written  by  Frederick  Sauer)  which 
produces  the  iteration  parameters  when  supplied  with  the  two  intervals  and  the  desired 
polynomial  degree. 


2.  The  extremal  polynomial  The  papers  mentioned  above  all  use  Chebyshev  polynomials 
in  some  essential  way,  so  we  first  note  that,  in  general,  the  required  polynomial  Pn  is 
unrelated  to  the  Chebyshev  polynomial  Tn  for  S  .  This  is  seen  in  the  analysis  of 
Atlestam  [1977]  or,  more  directly,  from  the  fact  shown  below  that  Pn  alternates  one  less 


To  recall,  the  Chebyshev  polynomial  Tn  for  the  compact  set  S  is  the  polynomial  of 

the  form  tn  +  En  \  B.t"'  which  is  as  small  as  possible  on  S  .  In  other  words,  T„  is  the 
1=0  1  n 

error  in  the  best  approximation  on  S  to  tn  from  {1^  ^0  th  .  By  contrast,  we  are 

interested  in  the  polynomial  Pn  which  is  the  error  in  the  best  approximation  on  S  to 

1  from  { B  .  t ^  } 

1  1 

We  now  reformulate  this  problem  as  follows.  Let  A  be  the  linear  functional  on  it 

n 

(  :=  the  polynomials  of  degree  n  or  less)  whose  value  at  p  is  p C 0 )  .  In  symbols, 

A  :  r  — -►  R  :  p  | — +  p(0 ) 
n 

An  extremal  for  A  is  any  polynomial  of  norm  1  at  which  A  takes  on  its  norm,  i.e., 

any  p  z  with  llpll^  =  1  and  Ap  =  BAH  •  Here 

II  pH  -  Unit  :=  max  |p(s)| 

S  scS 
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and 


If  now  p* 

II XII 


a.  a  i  0  iff  t .  <  0  < 
3  3+1  3 

is  an  extremal  for  X  ,  then  we  have 

*  n+1  *  n+1 
Xp  =  E  ct.p  (t.)  <  E  |a. ||p  (t.)l 

j=1  3  3  j=1  3  3 


t 


3+1 


< 


(  E  |a  .  I  ]  Up  li 

3=1  3 


II  X  II 


so  equality  must  hold  throughout  this  relationship.  In  particular, 

*  ★ 

sign(Oj)p  (t^)  =  lip  II  =  1  ,  all  j  . 

This  pins  down  p*  uniquely  once  we  know  all  the  tj's  ■  Explicitly, 

.  n+1 

p  =  Z  sign( i . (0 ) )  t . 

3=1  3  3 

Thus,  for  n  >  1  ,  p*  is  not  just  a  constant,  therefore  Theorem  2.15  of  Rivlin  [1974! 
shows  that  the  points  t^,  . ..,  tn+1  are  uniquely  determined. 

If  now  0  lies  to  one  side  of  [t^,tn+^]  ,  then  it  follows  that  p  alternates  or. 
t .j ,  ...,  tn+1  ,  hence  p*  is  necessarily  a  multiple  of  the  Chebyshev  polynomial  for 
SO[t1(tn+1]  .  Further,  |p*(t)|  >  1  for  t  not  in  [t1(tn+1l  .  We  conclude  that  in  tne 
case  of  particular  interest  to  us,  namely  when  0  is  in  the  convex  hull  of  S  ,  there  must 
be  some  k  for  which 


<  0  <  tk+1 

This  shows  our  assertion  at  the  beginning  of  this  section  that,  in  general,  Pn  need  only 
alternate  on  n  points  in  S  .  Further,  for  t^  <  t  <  tk+1, 

*  n+1 

p  (t)  =  Z  (sign( i , (0 ) ) ) i . ( t)  =  Z  (signt l , ( t) ) ) Z . ( t) 

j=1v  3  3  3  3 

=  E  U.  (t)  |  >  I  -  i  (t)  I  =  1 

3  3 


if  n  >  1  .  Consequently,  then 

t  =  b  :=  max  S  fl(-«,0J  ,  and  t,  -  c  :  =  min  S  f 0 , °°  1 
k  k+1 

We  gather  these  various  facts  in  the  following  theorem,  for  the  record. 


Wieorem  1  Assume  that  S  is  compact  and  does  not  contain  0  ,  and  n  >  1 .  Thor, 
;  a)  X  has  a  unique  canonical  representation,  and  =  n^  .tVtt^-t^l.  , 

j= 1 , . . . , n+1 . 


: 
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(b)  Correspondingly,  X  has  a  unique  extremal,  anti  this  is  given  by 


n+1  t  -  t. 

=  £  sign(f.(0))  1.  with  £.(t)  =  H  - -  ,  all  j. 

3=1  3  3  3  i*j  V  fci 


If,  in  addition,  0  is  in  the  convex  hull  of  S  ,  then 


(c)  tk  < 

0  for  some 

k  .  Further 

* 

,  j=1 , . ,k 

p  <tj) 

(-,3+1-k 

,  j=k+1 , . . • ,n+1 

and  p* ( t )  >  1 

for  tk  <  t  <  tk+1 

,  hence  * 

max  S  n and  tk+ ^  =  min  sn[(j,» 

We  pointed  out  earlier  that  Pn  =  p  /p  (0)  could  also  be  obtained  as  the  error  in  the 
Chebyshev  approximation  to  1  from  the  subspace  $.t^}  .  This  suhspace  forms  a  Haar 

set  on  S  as  long  as  S  does  not  contain  0  .  We  could  therefore  have  obtained  the  above 
characterization  from  general  criteria  such  as  Kolmogorov's  criterion,  but  the  derivation 
would  not  have  been  any  simpler.  We  note  that,  while  Pn  is  in  general  not  (a  multiple  of) 
the  Chebyshev  polynomial  for  S  ,  it  is  always  a  multiple  of  a  Zolotarev  polynomial  since 
its  alternations  over  n  points  characterize  it  as  the  error  in  a  best  approximation  to 


e  t"  +  e  ,t' 

n  n- 1 


n-1 


from  it 


n-2 


3.  Rentes  algorithm  for  the  extremal  polynomial  We  begin  with  a  statement  of  the 
algorithm.  In  it,  we  use  the  abbreviations 

b  :=  max  Sf~M-oo,0]  ,  c  :=  min  Sn  (0,M) 

introduced  earlier. 


Remes  algorithm  for  the  extremal  polynomial 

1.  Choose  t  *  (t^}j_|  in  S  ,  strictly  increasing,  and  with  tk  =  b,  +  ^ 

k  • 

2.  Set  p  :=  2"*Jsign(  (0  ) )  ,  with  £_.  ( t )  :=  ( t-t^  )  /( t . -t t)  /  all  j 

3.  Set  tfl  min  S  ,  tn+ 2  :=  max  S  and  construct  s  by 


c  for  some 


s 

D 


t.  for  j-k,k+1 
3 

the  first  of  the  possibly  two  maxima  of  p(t.)p  in 
for  j-1 , • . . ,k-1 ,k+2, • . • ,n+1 


‘Vi'Vi* 


s 


4.  Choose  t  from  s  as  follows: 

(a)  if  p(tQ)p(t1)  <  -1  ,  then  t  :=  (t0,  s^  ...,  sn),  and  increase  k  by  1  . 

( b)  if  P< bn+2 bn+i )  <  -1  ,  then  t  :=  (s1(  . ..,  sn+^,  tn+2 )  /  and  decrease  k  by  1 . 

( c)  otherwise,  t  :=  s 

5.  Set  t  :=  t  . 

6.  Iterate  steps  2  through  5  . 


Theorem  2  The  sequence  of  polynomials  produced  by  the  above  Remes  algorithm 

* 

converges  to  p  . 

Proof  We  first  note  that  the  algorithm  is  well  defined  at  step  4  in  that  only  one  of 
the  alternatives  (a)  and  (b)  is  possible.  Indeed,  if  both  (a)  and  (b)  were  to  occur,  then 
p  would  alternate  on  tQ,  ...,  tk,  tk+2  tn+2  ,  and  this  is  not  possible  for  a 

polynomial  of  degree  n  or  less  . 

As  to  the  convergence,  denote  by  p  the  polynomial  obtained  from  p  after  one 
iteration,  i.e.,  the  polynomial  constructed  from  the  sequence  t  obtained  at  step  4  .  We 
claim  that  1  <  p(t)  <■  p(t)  for  any  t  in  (b,c)  and  that  strict  inequality  holds  here 
unless  t  =  t  .  The  first  inequality  we  already  observed  earlier  (for  Pn).  As  to  the 
second,  we  have 


<-);’"k  p(t.) 

=  1  <  (-)j-kp(t.) 

for  j=1 , 2 , • . 

. ,  k 

-)  D(t  .  ) 

•  3 

j  -k  - 1  ~ 

=  1  <  (-)J  p(t.) 

for  j=k+ 1 ,  •  . 

•  ,  n+  1 

by  construction.  This  implies  that,  for  b  <  t  <  c  , 

~  n+1  ~  ~  ^ 

p(t)  -  p(t)  =  I  p(t.)  -  p(t.)lf.(t)  =  -  E lp(t  .  )-p( t  . ) I  I  4  .  ( t) I  <  0 
j  =  1  ’  3  13  3  3  1 

and  equality  occurs  only  if  p  =  p  on  the  n+1  points  t^,  ...,  t^+^  .  i.e.,  only  if 

P  ”  P  • 
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We  conclude  that  the  sequence  generated  by  the  Remes  algorithm  decreases  nonotonely  on 
(b ,c),  yet  is  bounded  below  there.  Hence  it  converges,  uniformly  on  any  finite  interval,  to 
some  polynomial  p*  .  This  polynomial  is  then  a  fixed  point  of  the  map  T  in  r ^  defined 
by 

T  :  p  - ►  p  , 

it  being  is  straightforward  to  show  that  T  is  continuous.  But  this  says  that  p*  is  the 
desired  extremal. 


4.  Efficient  computation  of  the  parameters  We  were  led  to  study  this  problem  by  the 
work  of  Roloff  [1979]  where  estimates  of  the  parameters  are  provided.  A  result  of  Roloff's 
states  that  the  zeros  of  Pn  are  approximately  distributed  in  S  in  a  proportion  which  is 
independent  of  n  .  One  might  hope  that  this  proportion  is  determined  by  measure,  i.e.,  a 
subinterval  of  S  containing  1/10  the  length  of  S  contains  about  1/10  of  the  zeros  of 
Pn  .  We  use  this  to  obtain  the  initial  guess  (in  step  1)  for  the  Remes  algorithm  but  we 
also  note  that  this  approximate  distribution  of  zeros  of  Pn  is  not  especially  good. 
Rather,  there  is  also  a  tendency  for  the  zeros  to  be  distributed  equally  among  the 
intervals  which  make  up  S  and  the  actual  distribution  resulting  from  these  conflicting 
tendencies  is  not  easily  predicted. 

The  Lagrange  basis  for  it  is  especially  suited  for  the  efficient  and  stable 

n 

implementation  of  the  Remes  algorithm  because  one  can  obtain  the  polynomial  p  associated 
with  the  current  point  sequence  t  without  any  computation,  because  the  basis  is  well 
conditioned  near  the  optimal  t  ,  and  because,  in  the  end,  the  zeros  of  Pn  =  p*/p*(0) 
are  particularly  easily  obtained  from  this  form. 

For  efficiency  in  evaluating  p  away  from  t  one  should  express  p  as 

n+1  n+1 

p(t)  «  It  (t-t  )  l  a./(t-t  ) 

J-1  j  j=1  3  j 

where 


a  :=  p( t .  )  /  II  (t  -t  ) 

j  ^  3  1 
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would  be  calculated  once  and  for  all.  Also,  the  derivative  of 


computed  by 


p*(t  )  «  H  (t  -t  )  Z  —r - * 

m  .  .  mi  t  -  t . 

1  *n  j  *m  m  j 

The  interior  local  extrema  of  p  are  estimated  by  parabolic  intern'  1  at i .  A- 
step,  the  unique  extremum  x*  ,  say,  of  the  parabola  matchino  p  at  t^,  t,  ar»  i  : 

found,  with  x  **  or  tj+^  depending  on  the  sian  of  p'(t^)  .  The  uni  ;ue  '-y-re~ 

* 

the  parabola  matching  p  at  tj,  tj,  and  x  is  then  taken  as  the  suitable-  a?.pr 
to  the  desired  local  extremum  of  p  .  This  is  a  version  of  the  st&ndar  :  techr.iaue  '<• 
locating  local  extrema  for  use  in  the  Remes  algorithm;  it  is  sufficiently  accor&i- 
quadratic  convergence  of  the  algorithm.  Note  that  in  our  particular  situation  th*. re- 
need  to  make  a  global  search  for  extrema  as  we  know  exactly  where  all  the  extrema  mi 
Once  the  extremal  polynomial  p*  is  found  sufficiently  accurately,  then  its  *-.• 
feund  by  the  Modified  Regula  Falsi.  The  zeros  are  already  bracketed  by  the  t-  ex.  v 
one  which  is  outside  the  interval  [t^,tn+1]  .  This  one  may  be  to  the  left  or  r *;«-'•  *■  . 

and  may  actually  be  at  infinity.  We  make  the  transformation  x  *  1  'v  f*:.4.  -* 
apply  the  same  method  to  [x^,xn4.^]  . 

Maehly's  second  method  (see  Maehly  [19631)  is  an  alternative  me  the  1  :*• 

Pn  .  Its  attraction  is  that  it  operates  directly  on  the  representation  or  ~ 
nn+1  , 


n“_'(1  -  z  t)  and  thus  does  not  need  the  second  phase,  the  computations  of  r' 


■judge  this  approach  to  be  less  efficient  overall  because  of  the  added  work  cf  s 
the  parameters  z j  of  the  new  polynomial  each  time  t  is  replaced  Vv  t  .  T*  ■ 
involves  the  solution  of  n+1  simultaneous  equations  with  a  full  coefficient  cat:: 
have  not  tried  this  approach;  see  Dunham  [1966]  for  some  remarks  concertina  irnr^v— 
convergence  of  this  method. 


5.  Properties  of  the  parameters  and  convergence  rate  of  the  iteration  An  t  . 

of  particular  examples  of  the  extremal  polynomials  shows  that  they  d^  not,  in  oor*  t 
belong  to  orthogonal  polynomial  families  and  that  they  are  n.y  ’  •  •  i  t  • 


w 


recurrence  relation.  It  follows  from  a  classical  result  of  Fekete  (see  Widom  [1969])  that 
HP  converges  as  n  tends  to  infinity  .  Consequently,  convergence  of  Richardson's 

first  order  method  with  these  parameters  is  geometric.  Of  course,  in  practice,  one  is  not 
likely  to  use  the  parameters  from  Pn,  PR+ ^ ,  Pn+2*  •••  *n  sequence,  but  is  likely  to  use 
the  parameters  for  a  fixed  n  cyclically.  One  still  obtains  geometric  convergence,  but 
examples  (such  as  seen  below)  show  that  n  should  be  rather  large  in  order  to  exploit  the 
method's  potential  fully. 

In  order  to  judge  the  convergence  rates  one  could  expect,  we  present  in  Figure  1  the 

graphs  of  p*  on  the  interval  '[b,c]  for  the  case  S  =  [-1,-.8]  [.2,1]  and  for  various 

an  ★ 

values  of  n  .  Recall  that  P  =  p  /p  (0)  ,  hence  IIP  II  =  1/p  (0)  .  Further,  it  is 

n  n 

worthwhile  at  this  point  to  realize  that  a  linear  change  in  the  independent  variable 

leaves  Pn  essentially  unchanged.  In  other  words,  if  this  particular  S  is  obtained  from 

some  interval  pair  S#  =  [a*,b*l  [c*,d*l  by  the  linear  change  t  «  f(t*),  so  that  -1  = 

f  ( a** ) ,  -.8  =  f(b*),  etc.,  then  the  polynomial  P*  for  S*  is  simply  P  -f  .  In 

n  n 

j*  *  * 

particular,  then  IIP  II  =  1/p  (f(0))  ,  thus  1/p  (t)  runs  through  the  possible  values  of 
n 

# 

such  IIP  II  as  t  runs  between  -.8  and  .2  .  Thus  as  one  moves  from  b  =  -.8  to 
n 

c  =  .2  ,  along  one  of  the  curves  for  fixed  n  ,  one  sees  the  effect  on  the  achievable  error 
reduction  of  the  location  of  the  origin  between  the  two  intervals  comprising  an  interval 
pair.  Note  that  the  rate  of  convergence  becomes  1  and  the  linear  system  becomes 
(possibly)  singular  as  f(0)  approaches  b  or  c  . 

Figure  2  shows  the  dependence  of  the  rate  of  convergence  on  the  relative  sizes  of  the 
two  intervals  which  make  up  S  .  A  contour  plot  is  given  of  the  maximum  possible  rate  of 
convergence  as  b  and  c  vary  while  a  =  -1  and  d  =  1  remain  fixed.  This  maximum  rate 
occurs  at  the  point  where  p*  is  at  a  maximum  (between  b  and  c)  which  depends  on  b 
and  c  .  This  rate  approaches  1  as  c-b  approaches  (1  and  becomes  ouite  fast  as  b 
and  c  approach  -1  and  -1,  respectively. 
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Figure  3  indicate*  the  effect  of  near  singularity  of  the  linear  system  on  the  rate  of 


convergence,  again  for  S  an  interval  pair  and  for  10  parameters,  we  plot  contour*  of  the 
logarithm  of  the  slope  of  p*  at  t  «  c  .  The  larger  this  slope,  the  less  the  rate  of 
convergence  is  degraded  by  the  origin  being  close  to  c  . 

Both  Figures  2  and  3  exhibit  a  somewhat  erratic  behavior  due  to  the  fact  that  the 
number  of  tj's  i-n  each  interval  is  a  discrete  function  of  h  and  c  .  This  suggests  that 
it  would  be  quite  difficult  to  obtain  accurate  and  simple  approximation  formulae  for  the 
parameter  distribution. 
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APPENDIX 


We  liat  hara  a  Fortran  program,  written  by  Frederick  Sauer,  which  realizes  the 
Rama ■  algorithm,  given  in  Section  3,  for  the  extremal  polynomial,  for  the  special 
case  that  S  consists  of  two  (nontrivial)  intervals.  A  more  general  version  which 
allows  for  S  to  consist  of  finitely  many  intervals,  some  or  all  of  which  may  even 
be  trivial,  has  also  been  written  by  Frederick  Sauer  who  will  submit  it  for 
publication  separately. 


PARAMETER  DEGREE-25 ,DEGP1«DEGREE+1 , DEGP3-DEGREE+3 
REAL  X( DEGP 3 ) , XT I LDE ( DEGP 3 ) , ALPH A( DEGP 1 ) , PROD ( DEGP t ) 

INTEGER  N, R , I PR I  NT , MAX IT 
DATA  I PRINT, MAX IT,  N  ,EPS2  ,EPS3  ,EPS4 
1  /  2  ,  20  ,  10.1.E-2  ,1.E-3  ,1.E-4/ 

CALL  SETUP (N,-1.,-.8,.2,1.,X,R) 

CALL  EXTREMt  X , N, R , I PR I  NT , MAX IT , EPS 2 , EPS 3 , EPS4 , XTILDE , ALPHA, 
1  PROD, PZERO) 

STOP 

END 


SUBROUTINE  EXTREM  { X , N, R , IPRINT , MAXIT , EPS2 , EPS3 ,EPS4 , XTILDE , 
1  ALPH A, PROD, PZERO) 


C  THIS  SUBROUTINE  FINDS  A  POLYNOMIAL  WHICH  IS  AN  EXTREMAL  OF  THE  LINEAR 

C  FUNCTIONAL  WHICH  IS  THE  POINT  EVALUATION  AT  THE  POINT  ZERO.  THE 

C  POLYNOMIAL  IS  CHOSEN  FROM  THE  SPACE  OF  POLYNOMIALS  OF  DEGREE  LESS 

C  THAN  OR  EQUAL  TO  N,  WHERE  N  IS  A  PARAMETER  SET  BY  THE  USER.  THE  NORM 

C  ON  THIS  SPACE  IS  DEFINED  TO  BE: 

C 

C  NORM(P)  -  MAX(P(T)  :  <X<1>  .LE.  T  .LE.  X(R+1))  OR 

C  <X(R+2)  .LE.  T  .LE.  X (N+3 ) ) ) 

C 

C  •••**  INPUT  •*•** 

C  X  -  AN  ARRAY  OF  DIMENSION  N+3  WHICH  CONTAINS  THE  STARTING  VALUES 
C  OF  THE  ABSCISSAE  IN  X(2>,  X ( 3  ) ,  ...  ,  X(N+2).  IN  X ( 1 )  AND 

C  X ( N+3 )  SHOULD  BE  THE  ENDPOINTS  (I.E.  A  AND  D  RESPECTIVELY) 

C 

C  N  -  THE  DEGREE  OF  THE  POLYNOMIAL. 

C 

C  R  SUCH  THAT  X(R+1)-B  AND  X(R+2)-C  IN  THE  PROVIDED  ARRAY  X. 

C  X< 1 ) ,X(R+1 ) ,X(R+2) ,X(N+3)  DEFINE  THE  INTERVALS  FOR  WHICH  THE 

C  NORM  OF  THE  POLYNOMIAL  IS  CALCULATED.  SEE  DEF.  OF  NORM  ABOVE. 

C 

C  IPRINT  -  --1  SUPPRESS  ALL  PRINTING.  PZERO  WILL  NOT  BE  CALCULATED 

C  -0  PZERO  -  P(0)  WILL  BE  CALCULATED  AND  PRINTED. 

C  -1  ALSO,  THE  ROOTS  OF  THE  POLYNOMIAL  WILL  BE  CALCULATED 

C  AND  PRINTED 

C  -2  ALSO,  FOR  EACH  ITERATION  A  MESSAGE  WILL  BE  PRINTED 

C 

C  MAXIT  -  THE  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED. 

C 

C  EPS 2  -  A  TOLERANCE  USED  TO  CONTROL  THE  MAJOR  ITERATION.  THE 

C  ROUTINE  WILL  STOP  WHEN  THE  MAXIMUM  CHANGE  OF  ANY 

C  PARTICULAR  ELEMENT  OF  X  CHANGES  BY  LESS  THAN  EPS2. 

C 
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EPS 3  ,  EPS4  -  TOLERANCE  PARAMETERS  USED  BY  SUBROUTINE  REGULA. 

THE  CALCULATED  VALUE  OF  THE  ROOT,  XRT,  WILL  BE 
SUCH  THAT  XRT  WILL  BE  WITHIN  EPS 3  OF  AN  ACTUAL 
ROOT  OR  ABS (P (XRT ) )  WILL  RE  LESS  THAN  EPS4. 


WORKSPACE  ***** 

XTILDE  -  AN  ARRAY  OF  DIMENSION  AT  LEAST  N+3. 
ALPHA  -  AN  ARRAY  OF  DIMENSION  AT  LEAST  N+1 . 
PROD  -  AN  ARRAY  OF  DIMENSION  AT  LEAST  N+1 . 


OUTPUT  ***** 

-  the  FINAL  VALUES  OF  THE  ABSCISSAE  OF  THE  EXTREME  POINTS 
CALCULATED  BY  THIS  ROUTINE. 

-  THE  FINAL  VALUE  OF  R  AS  DESCRIBED  ABOVE. 

ALPHA,  PROD  -  THE  FINAL  VALUES  OF  ALPHA  AND  PROD  AS  CALCULATED 
BY  SUBROUTINE  ALCALC  ARE  RETURNED  IN  THESE  ARRAYS  IF 
IPRINT  IS  NOT  EQUAL  TO  -1. 

PZERO  -  VALUE  OF  EXTREMAL  POLYNOMIAL  AT  0. 

XTILDE  -  CONTAINS  THE  RECIPROCALS  OF  THE  ROOTS  OF  THE  POLYNOMIAL  IF 
IPRINT-1  OR  2.  XTILDE ( 1  )  CORRESPONDS  TO  THE  ROOT  LYING  OUT¬ 
SIDE  THE  INTERVAL.  IT  WILL  BE  SET  TO  ZERO  IF  THE  ADDITIONAL 
ROOT  DOES  NOT  EXIST. 


*****  NOTE:  THIS  SUBROUTINE  EXPECTS  A  STARTING  GUESS  FOR  THE  VALUES 
OF  X,  AND  A  VALUE  FOR  R.  THIS  CAN  BE  ACCOMPLISHED  BY  A  CALL  TO 
SUBROUTINE  SETUP. 

INTEGER  R.N.NP1,  ITER,  IPRINT, MAXIT,  ISET,  NP3  ,NP2  ,  RP1  ,RP2 
REAL  X( 1 ) , XTILDE ( 1 ) , ALPHA)  1 ) ,PROD< 1 ) 

INTEGER  I,J 

REAL  EPS2, SIGN, SGSL, ERR, EPS3, PZERO 

NP1  -N+1 

NP2  -  N  +  2 

NP3-N+3 

ITER  -  0 

IF  (N  .LE.  1)  RETURN 

40  IF  (IPRINT  .EQ.  2)  WRITE  (6,640)  ITER, (X (I ) ,I-2,NP2 ) 

640  FORMAT( '0 AFTER  ',13,'  ITERATIONS  THE  POINTS  X  ARE* ,/ ( IX , 1 0E1 3 . 5 ) ) 
RP1-R+1 
RP2-R+2 

GET  THE  INTERPOLATING  POLYNOMIAL 

CALL  ALCALC(X,NP1,R, ALPHA, PROD, XTILDE) 

XTILDE ( 1 )— X( 1 ) 

SIGN- 1 . 

IF  ( MOD ( R , 2  )  .EQ.  1)  SIGN— 1 . 

XTILDE (NP 3 )-X (NP3 ) 

ISET  -  0 

IF  ISET  IS  CHANGED  THE  RESULTING  X'S  MUST  BE  SHIFTED. 

IF  (X < 1 )  .EQ.  X ( 2 ) )  GO  TO  45 

CALL  EVAL  (X ( 1  )  ,X,  ALPHA, NP1  ,VAL) 

IF  <VAL*SIGN  .GT.  1.)  ISET— 1 
45  DO  80  I-2.NP2 

IF  ((I  .EQ.  RP1 )  .OR.  (I.EQ.RP2))  GO  TO  70 
J-I-1 

CALL  OERIV(X, J, PROD, ALPHA, NP1 , SLOPE) 

SGSL— 1 • 
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IF  (SLOPE'SIGN  .GT.  0.)  SGSL— 1. 

IF  (SGSL  . EQ.  1)  U-I+1 
IF  (SLOPE  .  NE.  0.)  GO  TO  60 
XTILDE (I )-X(I ) 

GO  TO  80 

60  SLOPE— ABS(  SLOPE) 

CALL  INTERPd  , J.SIGN.X,  ALPHA,  W1  .SLOPE  ,SG  L ,  XTILDE  ’  I  1 

SIGN-SIGN 

GO  TO  80 

70  XTILDE ( I ) -X ( I ) 

SIGN-1. 

80  CONTINUE 

IF  (X(NP2 )  .EQ.  X ( NP3 ) )  GO  TO  85 

CALL  EVAL(X  (N+3),X,  ALPHA,  1 ,  VAI. ) 

IF  (VAL*SIGN  .LT.  -1.)  ISFT-1 
C  SET  UP  X  FOR  NEXT  ITERATION 
85  R  -  R  -  ISET 
ERR— 0 . 

DO  90  I  -  2,NP2 
J  =  ISET  +  I 

DIFF  -  ABS(  X(I)  -  XTILDE ( J )  ) 

IF  (DIFF  .GT.  ERR)  ERR  =  DIFF 
90  X ( I )  =  XTILDE ( J ) 

IF  (  ERR  .LT.  EPS 2 )  GO  TO  100 
ITER  ■  ITER  +  1 
IF  (ITER  .LT.  MAXIT)  GO  TO  40 
WRITE  (6,650)  MAXIT 

650  FORMAT( '0MAXIMUM  NUMBER  OF  ITERATIONS  EXCEEDED ,  “ACT 
100  IF  (IPRIOT  .EQ.  2)  WRITE  (6,660)  (X ( I ) , 1-2, NP2 ■ 

660  FORMAT! '0 AFTER  FINAL  ITERATION  THE  POINTS  X  APE', /(IX, 1 
IF  (IPRINT  .EQ.  -t)  RETURN 
C  EVALUATE  EXTREMAL  POLYNOMIAL  AT  0  . 

CALL  ALCALC( X , NP 1 , R , ALPHA , PROD , XTILDE ) 

CALL  EVAL(0.  ,X,  ALPHA,  NP1  ,  PZERO) 

WRITE ( 6,670 )  P7ER0 
670  FORMAT ( *  OP ( 0 )  =',E20.8> 

IF  (IPRINT  .EQ.  0)  RETURN 

C  FIND  THE  RECIPROCALS  OF  THE  ZEROS  OF  THE  EXTREMAL  POLY',. \ 
CALL  FINDZR(X, ALPHA, NP1,R,EPS3, EPS4 , XTILDE) 

RETURN 

END 


SUBROUTINE  ALC ALC ( X , M , R , ALPH A , PROD , WORK ) 

C  THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENTS  ALPHA  ANT  FRi'D 
C  USED  TO  EVALUATE  THE  LAGRANGE  INTERPOLATING  POLY NON  I  A!  »NP 
C  DERIVATIVES 


C 

C  *»»**  INPUT  ***** 

C  X  -  THE  ABSCISSAE  OF  THE  POINTS  TO  WHICH  THr  P  LY'I  >!':■•;  - 

C  FIT  APE  STORED  IN  X(2),  X ( 3  )  ,  X(4> . 

C  M  -  THE  ORDER  (DECREE  +  1  )  OF  THE  INTERPOL  AT  IN  ■  r'TV: 

C  R  -  IS  SUCH  THAT  X ( V* 1 )  =  P,  X(P+?>  =  C 

C 

C  ***••  OUTPUT  ***** 

C  ALPHA  -  AN  ARRAY  OP  DIMENSION  M  WHICH  IS  I'cp"  r  ; —  -  r  • 
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ALPHA(i) 


C 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 


M 

=  P(X<I+t))/  PRODUCT  (X(I+1)-X(J+1)) 

J»1 


J.NE.I 

WHERE  P(X)  TP.  THE  VALUE  OF  THE  POLYNOMIAL  AT  X.  IN  PARTICULAR, 

IF  I  .LE.  R 

PIXII  +  1  )  )- 

(-1 )*#(I+1-F)  IF  I  .CT.  R 

PROD  -  AN  ARRAY  OF  DIMENSION  M  WHICH  IS  USED  TO  STORE  THE  VALUES 

M 

PROD(T)  *  PRODUCT  (X(I+1)  -  X(J+1)) 

J-l 

J.NE.I 


****  W^RKSPACF  ***** 

WORK  -  AN  ARRAY  Oh  DIMENSION  M. 

INTEGER  R,M,1,J,MM1 

REAL  X  ( 1  )  , ALPH A(  M )  , PROD ( M ) , WORX( M) ,XI ,P, SIGN 
MM1  *  M- 1 
SIGN  -  1. 

IF(MOD( R , 2 )  . RQ.  0)  SIGN  -  -1. 

DO  10  1-1 ,  MMl 
10  WORK ( I )  -  X ( I +2  ) 

DO  30  J-t,M 
P  -  1. 

XI  -  X( J41  ) 

DO  20  1*1, MMl 

20  P  -  P* ( X I -WORK { I  )  ) 

PROD(J)  =  P 
WORK ( J )  *  XI 

IF  (J  .  FQ.  P+1)  SION  *  -  SIGN 
ALPHA(J)  -  STGN/P 
30  sign  -  -  srnu 
RETURN 
END 


SUBROUTINE  EVALfT,  v.  AT.PHA,M,P) 

C  THIS  SUBROUTINE  EVALUATES  THE  VALUE  OF  THE  LAGRANGE  INTERPOLATING 
C  POLYNOMIAL  AT  THE  POINT  T. 

C 

C  *****  INPUT  ***** 

C  T  -  POINT  AT  WHIrH  POLYNOMIAL  IS  TO  BE  EVALUATED. 

r  X  -  ARRAY  r*ONTA I  NT POINTS  AT  WHICH  THF  POLYNOMIAL  IS  KNOWN. 

r  AI.PHA  -  ARPAY  OF  DIMENSION  M  WHICH  HAS  THE  COEFFICIENTS  CALCULATED 
C  BY  SUBROUTINE  ALCALC. 

M  -  ONE  PI  up  THF  DFGRFF  OF  THE  POLYNOM I  A3.. 


^  •##*#  OUTPUT  ***** 

C  P  -  THE  VALUE  K  THE  POLYNOMIAL  A’l  T. 

r 

r 

INTEGER  Y,I 

PFM.  X  ( 1  )  ,  ALPPAf  M)  ,T.  r,C,DIFF,SGV 

'■  FVAI.ftVL  Tier  P^LYN'iMZAL 


-1- 


c 

C  =  1. 

P  =  0. 

DO  10  I  =  1,M 

DIFF  =  (T-X ( 1+1 ) ) 

IF  (DIFF  .EQ.  0.)  GO  TO  20 
C  »  C*DIFF 

10  P  =  P  +  ALPHA! I) /DIFF 
P  =  P*C 
RETURN 
20  SGN  *  1. 

IF  (MOD(M-I,2)  .EQ.  1)  SGN  =  -1. 
P  =  SIGN  (1.,  ALPHA(I>*SGN) 
RETURN 
END 


SUBROUTINE  OERIV(X , I , PROD , ALPHA, M , D ) 

C  THIS  SUBROUTINE  EVALUATES  THE  DERIVATIVE  OF  THE  LAGRANGE  INTER- 
C  POLATING  POLYNOMIAL  AT  THE  ITH  POINT  OF  X. 

C 

C  *****  input  ***** 

C  X  -  ARRAY  CONTAINING  POINTS  AT  WHICH  THE  POLYNOMIAL  IS  KNOWN. 

C  I  -  THE  DERIVATIVE  IS  TO  BE  EVALUATED  AT  X(I+1> 

C  PROD  -  AN  ARRAY  OF  DIMENSION  M  WHICH  CONTAINS  COEFFICIENTS  CALCU- 
C  LATED  BY  SUBROUTINE  ALCALC. 

C  ALPHA  -  AN  ARRAY  OF  DIMENSION  M  WHICH  CONTAINS  COEFFICIENTS 
C  CALCULATED  BY  SUBROUTINE  ALCALC. 

C  M  -  ONE  PLUS  THE  DEGREE  OF  THE  POLYNOMIAL. 

C 

C  *****  OUTPUT  ***** 

C  D  -  THE  VALUE  OF  THE  DERIVATIVE  AT  X(I+1). 


C 

INTEGER  I,M,J 

REAL  X(1),PROD(M),D,XS,ALPHA(M),AS 
XS  =  X(I+1 ) 

AS  =  ALPHA(I) 

D  =  0. 

DO  10  J=  1 ,  M 

IF  (J  .EQ.  I)  GO  TO  10 
D  =  D  +  (AS  +  ALPHA!  J)  )/(XS-X(  J+1  )  ) 
10  CONTINUE 
D  =  D*PROD(I) 

RETURN 

END 


SUBROUTINE  INTERP  (I ,J, SIGN, X, ALPHA, M, SLOPE, SGSL, XMIN) 

C  THIS  SUBROUTINE  FINDS  AN  APPROXIMATION  TO  XMIN  WHERE  XMIN  IS  SUCH 
C  THAT: 

C 

C  SIGN*P(XMIN)=MIN(SIGN*P(T):MIN(X(I) ,X(J)) .LE.T.LE.  M AX (X ( I ) ,X ( J) ) ) 

C 

C  WHERE  P (T )  IS  THE  INTERPOLATING  POLYNOMIAL  AT  THE  POINT  T. 

C 
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C  *****  INPUT  ***** 

C  I,J  -  INTEGERS  SO  THAT  X(I)  AND  X(J)  DEFINE  THE  ENDPOINTS  OF  THE 
C  INTERVAL  TO  BE  SEARCHED  FOR  XMIN. 

C  SIGN  -  REAL  VALUE,  EITHER  +1.  OR  -1.  SO  THAT  SIGN*P(X(I))  =  -1. 

C  X  -  ARRAY  OF  DIMENSION  M+2  WHICH  CONTAINS  THE  ABSCISSAE  OF  THE 
C  POINTS  USED  FOR  LAGRANGE  INTERPOLATION. 

C  ALPHA  -  COEFFICIENTS  FOR  LAGRANGE  POLYNOMIAL. 

C  M  -  ORDER  (DEGREE  PLUS  ONE)  OF  INTERPOLATING  POLYNOMIAL. 

C  SLOPE  -  -<ABS(P'(X(I)))) 

C  SGSL  -  HAS  THE  FOLLOWING  VALUE 

C  +1.  IF  X(I)  .LT.  X(J) 

C  -1.  IF  X ( I )  .GT.  X(J) 

C 

C  *****  OUTPUT  ***** 

C  XMIN  -  THE  APPROXIMATION  OF  XMIN  CALCULATED  BY  THE  ROUTINE 
C 

C  *****  METHOD  ***** 


C  LET  P(T )  BE  THE  LAGRANGE  INTERPOLATING  POLYNOMIAL  AT  THE  POINT  T. 

C  S,  THE  FIRST  APPROXIMATION  OF  XMIN,  IS  THE  MINIMUM  OF  THE  PARABOLA 

C  DETERMINED  BY  P(X(I)),  P<X(J)),  AND  P'(X(I)).  XMIN  IS  THEN  THE 

C  THE  POINT  AT  WHICH  THE  PARABOLA  INTERPOLATING  P(X(I)),  P'(X(D), 

C  AND  P(S)  TAXES  ON  ITS  MINIMUM. 

C 

INTEGER  I, J,M 

REAL  SIGN,X(1 ) ,  ALPHA(M) , SLOPE , SGSL , XMI N 
REAL  S,XMAX,F,XI ,XT,DIFF2 
XI  »  X (I ) 

IFIRST-1 

S*  ABS (XI  -  X(J) ) 

XMAX  =  S 
F-  1. 

C  CHECK  IF  END  POINT 

IF  ((J  .NE.  1)  .AND.  (J  .NE.  M+2))  GO  TO  10 
IF  (X(J)  .NE.  XI)  GO  TO  5 
XMIN-XI 
RETURN 

C  CALCULATE  INTERPOLATING  POLYNOMIAL  AT  X(J) 

S  XT-X(J) 

CALL  EVAL(XT,X, ALPHA, M,F) 

F-F*SIGN 

C  USE  FOR  FIRST  STEP  THE  SLOPE  AT  X(I)  AND  THE  POINT  S. 

10  DIFF2-  ((F+1.J/S  -  SLOPEJ/S 
IF  (DIFF2  .LE.  0)  GO  TO  999 
S*  -SLOPE/DIFF2/2. 

IF  (S  .LE.  XMAX)  GO  TO  30 
XMIN  -  SGSL*XMAX  +  XI 
RETURN 

30  XT*  SGSL*S  +  XI 

CALL  EVAL(XT,X,  ALPHA, M,F) 

F*  SIGN*F 

IF  (I FIRST  .NE.  1)  GO  TO  999 
IFIRST-2 
XMAX  -  S 
GO  TO  10 

999  XMIN-  XI  +  SGSL*S 
RETURN 
END 


SUBROUTINE  FINDZP ( X , ALPH A , M , P , EPS 3 , EPS4 , WORK ) 

C  THIS  SUBROUTINE  FINDS  THE  ZEROS  OF  THE  LAGRANGE  INTERPOLATING  POLV- 

C  NON  I AL  USING  THE  MODIFIED  REGULA  FALSI  METHOD.  THE  RESULTS  ART' 

C  PRI  TED  OUT. 
c  *****  INPUT  ***** 

C  X  -  ARRAY  CONTAINING  THE  ABSCISSAE  OF  THE  DATA  POINTS. 

C  ALPHA  -  ARRAY  CONTAINING  COEFFICIENTS  CALCULATED  BY  SUBROUTINE 

C  ALCALC. 

C  M  -  ORDER  OF  LAGRANGE  POLYNOMIAL. 

C  R  -  INTEGER  SUCH  THAT  ZERO  LIES  IN  THE  INTERVAL  X(P+1)  TO  X*P»2.. 

C  EPS3.EPS4  -  TOLERANCE  PARAMETERS  USED  BY  SUBROUTINE  REGULA.  "UK 
C  CALCULATED  VALUE  OF  THE  ROOT,  XPT,  WILL  EE  SUCH  THAT 

C  XRT  WILL  BE  WITHIN  EPS  3  OF  AN  ACTUAL  ROOT  OR 

C  ABS ( P  t XRT) )  WILL  BE  LESS  THAN  EPS4. 

C 

C  *****  WORKSPACE  ***** 

C  WORK  -  AN  ARRAY  OF  DIMENSION  AT  LEAST  M. 

C 

INTEGER  M,R,J,I,MM1 

REAL  X( 1 ) ,WORX(M) , ALPH A( M) , SIGN, EPS3 , EPS4 

REAL  A, B , F A, FB , SIGN1 , SUM, XRT 

EXTERNAL  EVAL , EVALI 

MM1 -M- 1 

SIGN  =  1 . 

IF  f  MOD ( R , 2 )  .EQ.  1)  SIGN--1 • 

SIGN1  =  SIGN 
DO  20  J=2,M 

IF  (J  .EQ.  R+1 )  GO  TO  20 
A  -  X(J ) 

B  =  X( J+1 ) 

FA  =  -SIGN 
FB  =  SIGN 

CALL  REGULA( EVAL, A,B ,F A,FB ,X , ALPHA,M , EPS 3 , EPS4 , XRT ) 

SIGN=-SIGN 
I  =  J 

IF  {J  .GT.  R+1)  I=J-1 
WORK ( I )  =  XRT 

20  CONTINUE 

WRITE  (A, 620)  f  WORK.  { I  )  ,  I  ~2  ,MM  1  ) 

620  FORJMATCOTHE  ROOTS  LYING  IN  THE  INTERVAL  (X  (  2  )  ,  X  I M+ 1  )  )  ARE', 

1  /(  1X,5F.20.7)) 

DO  21  J=2,M 

21  WORK(J)  =  1 • /WORK ( J ) 

C  FIND  ROOT  OUTSIDE  OF  INTERVAL  IF  IT  EXISTS 
SUM  =  0. 

ALMAX  =  0. 

DO  30  1=1, M 

ARSAL  =  ABS{ ALPH A( I ) ) 

IF  (ABSAL  .GT.  ALMAX)  ALMAX  *  ABSAL 
30  SUM  =  SUM  +  ALPH A(I) 

IF  1 ABS (SUM!  .GT.  ALMAX  *  1.E-4)  GO  TO  40 
36  WRITE ( 6,636 ) 

636  FORMAT* ’OTHE  ADDITIONAL  ROOT  CANNOT  BE  ACCURATELY  O  \i  AT'"  ' 

1  '  OR  IS  INFINITE 1  ) 


m 


r — 


Wi  •!  - (  1  C  . 

RETURN 

40  tk  (sum*lugn  .0*!.  ;).)  G'*  TO  70 
C  THE  ADDITIONAL  ROOT  IS  TO  THE  LEFT  OF  X ( 2 ) 

A  -  1 • /X ( 2  ) 
n  =  - i . e-4 
FA  -  -SICN1 

CALL  EVALI  ( B , X , ALPHA , M , FB ) 

IF  (FA*FB  .GE.  r . )  GO  TO  35 
GO  TO  90 

C  THE  ADDITIONAL  IS  TO  THE  RICHT  OF  X(MM  ) 

70  P  -  1 . /X ( M+ 1 ) 

A  =  1 .E-4 
F3“~SICN 

CALL  EVAi.l  l  A#  a,  tl ,M , FA ) 

IF  { F .V FB  .GE.  0  )  G-'<  TO  >r> 

90  CALL  REGUL.M  EVAL1 ,  A,  L,  rA,  FB ,  X,  ALPHA,  M,  EPS3,EPS4,  XLT) 
WRITE  (6,097)  XPT 
WORK  ( 1  )sXRT 

69C  FORMATCOTHE  RECIPROCAL  OF  THE  ADDITIONAL  ROOT  TS  A "  2". 

RETURN 
END 


S  UP  R  OUT  INF  F.EGtJL  A  {  FCT ,  A,  H ,  FA,  FB ,  X  ,  ALPH  A ,  M  ,  EPS 1 ,  FPS4  ,  XRT ) 

C  THIS  CT.TRP.0UT INC  FINDS  THE  ZERO  OF  THE  LAGRANGE  POLYNOMIAL  WHICH 
C  IN  THE  INTERVAL  (A,P>.  THE  MODIFIED  RECTI  LA  FALSI  METHOD  IS  USED. 

C  *****  INPUT  ***** 

C  A  -  LEFT  END  POINT  OF  INTERVAL  TN  WHICH  A  ROOT  LIES. 

C  D  -  RIGHT  END  POINT  CF  INTERVAL  IN  WHICH  A  ROOT  LIES. 

C  FA  -  THE  VALUE  OF  THE  POLYNOMIAL  AT  THE  POINT  A. 

C  FB  -  THE  V ALU*..'  OT  ’HYNOMIAL  AT  THE  POINT  B. 

C  X  -  ARRAY  CONrATMTNO  THE  ABSCISSAE  OF  THE  DATA.  POINTS. 

C  ALPHA  -  ARRAY  CONTAINING  COEFFICIENTS  CALCULATED  BY  SUBROUTINE 

C  AT. 

C  M  -  OHDEH  CT  LAOL'.NGf  , ■■ i  "'PV  J .  AL  . 

C  EPS 3 , EPS4  -  TOLERANCE  PARAMETERS.  THE  CALCULATED  VALUE  OF  XRT 
C  WILL  Bfc  SUCH  1  HAT  XRT  WILL  UE  WITHIN  EPS >  OK  AN 

C  ACTUAT.  FOOT  OF.  ARS  ( P  ( XF.T  )  }  WILL  BE  LESS  THAN  EPS 4  . 

r 

q  *****  output  ***** 

C  XRT  -  THE  CAI.OUI.ATFD  VALUE  OF  THE  FOOT. 

C 

INTEGER  '* 

REAL  X(1  )  ,  ALPHATM)  ,  A,R , p  A,  FB ,  EPS3  ,  t'PS4 , 1-V ,  XPT ,  FXFT 
EXTERNAL  FCT 
F’W  -  FA 

10  XRT  '  (op* A  -  F A*n)/(^P  ~  FA) 

OALT.  FCT  (  XRT,  X  , ALPHA, M ,  I-  XRT  ) 

IF  (ARS(FXPT)  .LE.  EPS4 )  RETURN 
TF  (  S T ON ( 1 . , FXPT ) * F A  .OT.  0.  )  GO  TO  12 

H  =  XRT 
FP  -  FXF'" 

IF  (S  IGN(  1  .  ,  FXPT  )  *  F»v  .GT.  0.)  FA-FA/2  • 

GO  To  15 
12  A  -  XRT 
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FA  =  FXRT 

IF  (SIGN( 1 . ,FXRT)*FW  ,GT*  0.)  FB=FB/2 . 
15  FW  =  FXRT 

IF  (B-A  .GT.  EPS 3 )  GO  TO  10 
XRT=  (A  +  B)  /  2. 

RETURN 

END 


SUBROUTINE  EVALI (T , X , ALPHA , M , P ) 

C  THIS  SUBROUTINE  EVALUATES  (T** ( M-1 ) ) *P ( 1 /T)  WHERE  P ( 1 /T )  IS  THE 
C  LAGRANGE  INTERPOLATING  POLYNOMIAL  EVALUATED  AT  1 /T 
C 

p  *****  INPUT  ***** 

C  T  -  POINT  AT  WHICH  POLYNOMIAL  IS  TO  BE  EVALUATED . 

C  X  -  ARRAY  CONTAINING  POINTS  AT  WHICH  THE  POLYNOMIAL  IS  KNOWN. 

C  ALPHA  -  ARRAY  OF  DIMENSION  M  WHICH  HAS  THE  COEFFICIENTS  CALCULATED 
C  BY  SUBROUTINE  ALCALC. 

C  M  -  ONE  PLUS  THE  DEGREE  OF  THE  POLYNOMIAL. 

C 

c  *****  OUTPUT  ***** 

C  P  -  THE  VALUE  OF  THE  POLYNOMIAL  AT  T. 

C 

c 

INTEGER  M, I 

REAL  X  ( 1  )  ,  ALPHA(M) ,T,P,C,DIFF,SGN 
C 

C  EVALUATE  THE  POLYNOMIAL 
C 

C  =  1. 

P  =  0. 

DO  10  I  »  1,M 

DIFF  =  ( 1 .-T*X(I+1 ) ) 

IF  (DIFF  .EQ.  0.)  GO  TO  20 
C  =  C*DIFF 

10  P  =  P  +  ALPHA(I )/DIFF 
P  =  P*C 
RETURN 
20  SGN  =  1. 

IF  (MOD( M-1 , 2 )  . EQ.  1)  SGN  =  -1  . 

P  =  SIGN  (1.,  ALPHA(I)*SGN)*T**(M-1 ) 

RETURN 

END 


SUBROUTINE  SFTUP  (N,  A,  B,  C,  D,  X,  R) 

C  THIS  SUBROUTINE  SETS  UP  THE  ARRAY  X  IN  PREPARATION  FOP  A  CALL  TO 
C  SUBROUTINE  EXT REM. 

C 

o  *****  INPUT  ***** 

C  N  -  THE  DEGREE  OF  THE  POLYNOMIAL  TO  RF  USED. 

C  A,  B ,  C,  D  -  INTERVAL  ENDPOINTS.  THE  INTERVAL  PAIR  IS  (A,B),  (C#D) 

C 

r  *****  output  ***** 

c  X  -  AN  APR  AY  OF  DIMENSION  N+3  CONTAINING  THE  CALCULATE!)  VALUES. 


o  o 


R  -  INTEGER  SUCH  THAT  X(R+1 )  -  B  AND  X(R+2)  -  C. 

INTEGER  N,R,RM1,NMR 

REAL  X(1),  A,B,C,D,DEL1 , DEL2 

IF  (N  .LE.  1)  RETURN 

DELI  =  B  -  A 

DEL2  -  D  -  C 

R  -  INT(  FLOAT (N-1 )  •  (DELI /(DELI  +  DEL2 ) ) )  +  1 
IF  (R  .LE.  1)  GO  TO  20 
RM1  “  R  -  1 

DELI  »  DELI  /  FLOAT ( RM 1 ) 

DO  10  1=2,  R 

10  X(I)  «  A  +  FLOAT(I-2)*DEL1 

20  NMR  =  N-R 

IF  (NMR  .LT.  1)  GO  TO  40 

DEL2  ”  DEL 2  /  FLOAT  (NMR) 

DO  30  1  =  1, NMR 

30  X(R+I+1)  =  C  +  FLOAT(I-1 )*DEL2 

40  X(1)  =  A 
X(R+1 )  =  B 
X(N+2)  =  D 
X ( N+3 )  =  D 
RETURN 
END 
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